function group_matrix =  penalized_hierarchical_clustering...
    (dissimilarity,G_vec,Ng_lb)
% This function takes a dissimilarity matrix and a vector of group numbers,
% and returns a sequence of clustering structures, using hierarchical
% clustering algorithm (average linkage) that penalizes imbalance among 
% group sizes. 

[n,~] = size(dissimilarity);
mvec = tril(true(size(dissimilarity)),-1); % lower triangular parts 
Rd = dissimilarity(mvec); % tranfer dissimilarity matrix into a vector 

% choose tuning parameter such that the smallest group is large enough
% b = .6/G_vec(end);
% Ng_lb = b*n; % smallest group size as a fixed proportion of n
lambda = 0;
step = .01;
hcRd = linkage_penalty(Rd','co',lambda);
group = cluster(hcRd,'maxclust',G_vec(end));
freq  = sum(repmat(unique(group),1,length(group))==repmat...
    (group',length(unique(group)),1),2);
Ng_min = min(freq);
while (Ng_min<Ng_lb)
    lambda = lambda+step;
    hcRd = linkage_penalty(Rd','co',lambda);
    group = cluster(hcRd,'maxclust',G_vec(end));
    freq  = sum(repmat(unique(group),1,length(group))==repmat...
        (group',length(unique(group)),1),2);
    Ng_min = min(freq);
    
    % break rule
    if lambda >= 1
        break;
    end
end

% lambda

% use the chosen tuning paprameter lambda
group_matrix = zeros(n,numel(G_vec));
for jj = 1 : numel(G_vec)
    group_matrix(:,jj) = cluster(hcRd,'maxclust',G_vec(jj));
end

